stanEx3
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}non linear regression with one explanatory variable
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}ex8-0.stan
single linear regression
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -5186.51
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -43.7983 0.0117007 0.00153088 1 1 40
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -43.80
## b0 3.31
## b1 2.13
## s 5.42
## m0[1] 11.66
## m0[2] 22.90
## m0[3] 42.49
## m0[4] 11.58
## m0[5] 19.80
## m0[6] 39.50
## m0[7] 26.65
## m0[8] 13.29
## m0[9] 21.96
## m0[10] 9.68
## m0[11] 24.23
## m0[12] 31.88
## m0[13] 26.46
## m0[14] 40.10
## m0[15] 33.99
## m0[16] 4.49
## m0[17] 40.67
## m0[18] 19.25
## m0[19] 6.71
## m0[20] 18.32
## y0[1] 14.41
## y0[2] 14.42
## y0[3] 44.48
## y0[4] 4.80
## y0[5] 20.06
## y0[6] 41.77
## y0[7] 31.97
## y0[8] 3.64
## y0[9] 17.72
## y0[10] 13.58
## y0[11] 21.38
## y0[12] 32.26
## y0[13] 21.65
## y0[14] 35.87
## y0[15] 42.02
## y0[16] -7.79
## y0[17] 40.27
## y0[18] 17.12
## y0[19] 3.28
## y0[20] 19.37
## m1[1] 4.49
## m1[2] 8.72
## m1[3] 12.94
## m1[4] 17.16
## m1[5] 21.38
## m1[6] 25.60
## m1[7] 29.83
## m1[8] 34.05
## m1[9] 38.27
## m1[10] 42.49
## y1[1] -4.85
## y1[2] 6.70
## y1[3] 11.99
## y1[4] 17.65
## y1[5] 29.19
## y1[6] 24.13
## y1[7] 20.32
## y1[8] 39.09
## y1[9] 39.25
## y1[10] 45.54
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -43.71 -43.36 1.36 1.07 -46.44 -42.27 1.01 563 555
## b0 3.30 3.18 2.74 2.66 -0.73 7.98 1.02 283 325
## b1 2.14 2.14 0.25 0.24 1.71 2.54 1.01 367 539
## s 6.18 6.04 1.15 1.03 4.63 8.35 1.00 1278 1122
## m0[1] 11.66 11.54 1.96 1.86 8.75 14.93 1.02 299 431
## m0[2] 22.91 22.89 1.39 1.28 20.61 25.26 1.01 788 1047
## m0[3] 42.53 42.57 2.68 2.57 38.08 46.80 1.00 859 1286
## m0[4] 11.58 11.47 1.96 1.88 8.66 14.87 1.02 299 431
## m0[5] 19.81 19.81 1.45 1.34 17.47 22.28 1.01 496 654
## m0[6] 39.54 39.60 2.39 2.29 35.60 43.32 1.00 1052 1234
## m0[7] 26.68 26.67 1.45 1.39 24.26 29.09 1.00 1765 1321
## m0[8] 13.30 13.19 1.82 1.72 10.53 16.40 1.02 310 435
## m0[9] 21.97 21.95 1.40 1.31 19.70 24.32 1.01 669 960
## m0[10] 9.68 9.57 2.13 2.01 6.50 13.29 1.02 292 401
## m0[11] 24.25 24.26 1.40 1.30 21.90 26.55 1.01 1065 1203
## m0[12] 31.91 31.94 1.73 1.66 29.01 34.81 1.00 2098 1391
## m0[13] 26.48 26.47 1.45 1.39 24.08 28.87 1.00 1721 1338
## m0[14] 40.14 40.19 2.44 2.35 36.07 44.02 1.00 1000 1250
## m0[15] 34.03 34.06 1.89 1.82 30.90 37.18 1.00 1875 1358
## m0[16] 4.48 4.37 2.62 2.55 0.61 8.94 1.02 284 325
## m0[17] 40.71 40.75 2.50 2.38 36.57 44.69 1.00 956 1284
## m0[18] 19.26 19.26 1.47 1.35 16.89 21.75 1.01 467 627
## m0[19] 6.71 6.60 2.40 2.32 3.15 10.75 1.02 286 357
## m0[20] 18.33 18.33 1.51 1.36 15.93 20.92 1.02 426 568
## y0[1] 11.44 11.42 6.54 6.34 0.90 22.24 1.00 1574 1696
## y0[2] 22.75 22.71 6.39 6.28 11.82 32.46 1.00 2067 1917
## y0[3] 42.60 42.66 6.70 6.63 31.83 53.94 1.00 1990 1839
## y0[4] 11.61 11.56 6.54 6.51 0.57 22.28 1.00 1250 1616
## y0[5] 20.05 20.00 6.53 6.29 9.09 30.64 1.00 2005 1981
## y0[6] 39.51 39.62 6.80 6.64 28.57 49.98 1.00 1896 1785
## y0[7] 26.71 26.83 6.48 6.09 16.03 37.37 1.00 1811 2013
## y0[8] 13.50 13.34 6.42 6.22 3.32 23.64 1.00 1809 1724
## y0[9] 22.14 22.11 6.42 6.20 11.51 32.87 1.00 1789 1930
## y0[10] 9.97 9.97 6.65 6.26 -0.64 21.04 1.00 1286 1727
## y0[11] 24.55 24.54 6.40 5.99 13.83 35.06 1.00 2118 1735
## y0[12] 32.04 32.07 6.70 6.52 21.19 43.09 1.00 2099 2025
## y0[13] 26.49 26.50 6.55 6.23 15.52 37.34 1.00 2088 1922
## y0[14] 40.34 40.40 6.93 6.54 28.95 51.26 1.00 1701 1914
## y0[15] 33.87 33.90 6.69 6.44 22.78 44.80 1.00 1893 1902
## y0[16] 4.56 4.43 6.96 6.54 -6.94 16.39 1.00 1032 1329
## y0[17] 40.75 40.62 6.98 6.55 29.26 52.18 1.00 1916 1681
## y0[18] 19.16 19.12 6.32 6.27 8.78 29.44 1.00 1818 1744
## y0[19] 6.78 6.76 6.47 6.36 -3.50 17.41 1.00 1267 1644
## y0[20] 18.30 18.19 6.68 6.53 7.77 29.45 1.00 1517 1945
## m1[1] 4.48 4.37 2.62 2.55 0.61 8.94 1.02 284 325
## m1[2] 8.71 8.60 2.21 2.09 5.42 12.43 1.02 289 384
## m1[3] 12.94 12.83 1.85 1.75 10.16 16.05 1.02 307 435
## m1[4] 17.17 17.14 1.57 1.43 14.68 19.86 1.02 383 408
## m1[5] 21.40 21.38 1.41 1.32 19.10 23.77 1.01 613 863
## m1[6] 25.62 25.60 1.42 1.34 23.22 27.97 1.00 1490 1275
## m1[7] 29.85 29.88 1.60 1.51 27.15 32.50 1.00 2164 1406
## m1[8] 34.08 34.11 1.90 1.82 30.94 37.25 1.00 1866 1356
## m1[9] 38.31 38.36 2.27 2.19 34.60 41.95 1.00 1179 1285
## m1[10] 42.53 42.57 2.68 2.57 38.08 46.80 1.00 859 1286
## y1[1] 4.58 4.59 6.94 6.61 -6.93 16.18 1.00 1101 1436
## y1[2] 8.54 8.59 6.68 6.45 -2.33 19.49 1.00 1457 1701
## y1[3] 12.67 12.50 6.39 5.90 2.17 23.60 1.00 1350 1587
## y1[4] 16.93 17.03 6.59 6.43 6.17 27.76 1.00 1645 1790
## y1[5] 21.53 21.58 6.29 6.04 11.10 31.80 1.00 1831 1918
## y1[6] 25.83 25.86 6.39 6.03 15.49 35.93 1.00 2147 1991
## y1[7] 29.65 29.71 6.60 6.30 19.23 40.12 1.00 2000 1817
## y1[8] 34.07 33.95 6.53 6.19 23.09 44.91 1.00 1850 1825
## y1[9] 38.08 37.96 6.79 6.35 27.22 49.30 1.00 1672 1851
## y1[10] 42.51 42.50 6.86 6.63 31.28 53.59 1.00 1831 1649
quadratic regression
y=b0+b2(x-b1)**2
ex8-1.stan
quadratic regression
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real b2;
real<lower=0> s;
}
model{
y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b2*(x[i]-b1)^2;
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b2*(x1[i]-b1)^2;
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)## Initial log joint probability = -89.708
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 28 -4.82174 0.000298523 7.80957e-06 1 1 39
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -4.82
## b0 4.63
## b1 4.04
## b2 0.56
## s 0.77
## m0[1] 5.18
## m0[2] 8.11
## m0[3] 5.53
## m0[4] 8.86
## m0[5] 4.63
## m0[6] 9.80
## m0[7] 5.15
## m0[8] 12.27
## m0[9] 7.66
## m0[10] 15.36
## m0[11] 6.99
## m0[12] 12.21
## m0[13] 4.88
## m0[14] 4.95
## m0[15] 8.12
## m0[16] 7.16
## m0[17] 4.94
## m0[18] 7.08
## m0[19] 10.05
## m0[20] 11.07
## y0[1] 6.13
## y0[2] 7.25
## y0[3] 6.19
## y0[4] 8.25
## y0[5] 5.36
## y0[6] 9.58
## y0[7] 6.34
## y0[8] 11.21
## y0[9] 6.26
## y0[10] 14.09
## y0[11] 6.99
## y0[12] 12.24
## y0[13] 4.74
## y0[14] 6.66
## y0[15] 8.64
## y0[16] 8.00
## y0[17] 5.07
## y0[18] 7.06
## y0[19] 10.41
## y0[20] 11.04
## m1[1] 12.21
## m1[2] 8.97
## m1[3] 6.62
## m1[4] 5.18
## m1[5] 4.63
## m1[6] 4.98
## m1[7] 6.23
## m1[8] 8.37
## m1[9] 11.42
## m1[10] 15.36
## y1[1] 12.68
## y1[2] 8.86
## y1[3] 7.05
## y1[4] 4.56
## y1[5] 3.89
## y1[6] 5.40
## y1[7] 6.13
## y1[8] 8.43
## y1[9] 10.26
## y1[10] 14.49
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -9.25 -7.12 6.83 1.65 -32.40 -5.49 1.11 24 11
## b0 4.65 4.64 0.49 0.30 4.08 5.18 1.03 309 289
## b1 2.76 4.04 4.86 0.08 -15.61 4.18 1.09 29 11
## b2 0.52 0.55 0.14 0.04 0.01 0.62 1.10 27 12
## s 1.14 0.90 0.86 0.19 0.68 3.76 1.12 37 11
## m0[1] 5.40 5.21 0.86 0.32 4.71 7.64 1.10 28 11
## m0[2] 8.16 8.13 0.47 0.32 7.56 8.78 1.05 1028 24
## m0[3] 5.73 5.57 0.81 0.31 5.07 7.73 1.11 28 11
## m0[4] 8.73 8.83 0.59 0.33 7.47 9.35 1.08 34 11
## m0[5] 4.87 4.67 0.92 0.31 4.12 7.32 1.11 28 11
## m0[6] 9.59 9.76 0.82 0.37 7.48 10.33 1.09 28 11
## m0[7] 5.37 5.19 0.87 0.32 4.68 7.63 1.11 28 11
## m0[8] 12.06 12.24 0.99 0.52 9.90 13.07 1.08 34 11
## m0[9] 7.74 7.68 0.50 0.31 7.16 8.39 1.05 87 16
## m0[10] 14.94 15.28 1.69 0.71 10.26 16.43 1.09 29 11
## m0[11] 7.01 7.00 0.37 0.27 6.50 7.51 1.05 697 58
## m0[12] 11.81 12.14 1.45 0.51 7.27 12.94 1.10 26 11
## m0[13] 5.12 4.93 0.91 0.32 4.40 7.51 1.11 28 11
## m0[14] 5.14 4.98 0.78 0.29 4.47 7.11 1.12 30 11
## m0[15] 8.04 8.10 0.44 0.30 7.32 8.59 1.06 67 16
## m0[16] 7.27 7.20 0.55 0.30 6.69 8.07 1.07 45 12
## m0[17] 5.14 4.98 0.79 0.29 4.47 7.11 1.12 30 11
## m0[18] 7.20 7.11 0.56 0.31 6.61 8.02 1.07 42 12
## m0[19] 9.83 10.01 0.88 0.38 7.46 10.61 1.09 27 11
## m0[20] 10.77 11.02 1.14 0.44 7.37 11.70 1.10 26 11
## y0[1] 5.45 5.24 1.67 0.98 3.57 7.76 1.04 77 15
## y0[2] 8.16 8.15 1.41 1.04 6.31 9.99 1.05 2090 61
## y0[3] 5.70 5.59 1.52 0.99 3.79 7.86 1.05 500 21
## y0[4] 8.72 8.84 1.65 1.03 6.61 10.58 1.06 972 26
## y0[5] 4.92 4.68 1.78 1.05 2.97 7.21 1.06 64 13
## y0[6] 9.58 9.73 1.64 1.03 7.54 11.49 1.05 172 19
## y0[7] 5.38 5.22 1.76 0.98 3.44 7.69 1.06 207 15
## y0[8] 12.08 12.24 1.69 1.07 10.01 13.94 1.05 459 26
## y0[9] 7.72 7.67 1.37 0.97 5.97 9.51 1.04 1888 102
## y0[10] 14.90 15.23 2.21 1.20 11.71 17.18 1.07 40 12
## y0[11] 7.05 7.00 1.45 0.97 5.24 8.99 1.04 1921 41
## y0[12] 11.78 12.10 2.14 1.13 9.08 13.77 1.07 58 12
## y0[13] 5.11 4.96 1.66 1.03 3.16 7.30 1.06 299 16
## y0[14] 5.13 5.00 1.56 0.99 3.27 7.42 1.06 255 21
## y0[15] 8.04 8.09 1.49 0.98 6.13 9.87 1.05 1898 52
## y0[16] 7.26 7.17 1.49 0.97 5.44 9.25 1.04 1758 35
## y0[17] 5.13 4.99 1.62 0.96 3.27 7.03 1.04 734 21
## y0[18] 7.16 7.14 1.48 0.99 5.26 9.00 1.05 1968 47
## y0[19] 9.84 9.99 1.68 1.03 7.74 11.81 1.06 248 18
## y0[20] 10.80 10.98 1.74 1.09 8.46 12.78 1.06 113 16
## m1[1] 11.81 12.14 1.45 0.51 7.27 12.94 1.10 26 11
## m1[2] 8.83 8.93 0.61 0.33 7.45 9.46 1.08 33 11
## m1[3] 6.67 6.64 0.40 0.27 6.16 7.30 1.06 244 19
## m1[4] 5.35 5.21 0.72 0.29 4.70 7.04 1.11 31 11
## m1[5] 4.87 4.67 0.91 0.31 4.12 7.28 1.12 29 11
## m1[6] 5.21 5.02 0.89 0.32 4.51 7.54 1.11 28 11
## m1[7] 6.40 6.27 0.69 0.30 5.78 7.79 1.10 31 11
## m1[8] 8.41 8.38 0.46 0.33 7.81 9.00 1.04 1504 67
## m1[9] 11.26 11.39 0.82 0.47 9.75 12.15 1.07 41 11
## m1[10] 14.94 15.28 1.69 0.71 10.26 16.43 1.09 29 11
## y1[1] 11.81 12.10 2.01 1.08 9.02 13.89 1.05 50 12
## y1[2] 8.85 8.91 1.52 1.00 6.89 10.71 1.06 1718 27
## y1[3] 6.70 6.66 1.51 1.00 4.78 8.60 1.05 1346 49
## y1[4] 5.38 5.24 1.58 1.00 3.53 7.52 1.05 569 20
## y1[5] 4.86 4.71 1.69 0.99 2.96 6.81 1.05 167 21
## y1[6] 5.24 5.05 1.67 1.00 3.32 7.46 1.05 115 21
## y1[7] 6.35 6.26 1.69 1.00 4.53 8.28 1.06 976 24
## y1[8] 8.46 8.40 1.52 0.99 6.64 10.37 1.05 2098 38
## y1[9] 11.26 11.39 1.72 1.06 9.15 13.27 1.05 552 27
## y1[10] 14.91 15.23 2.21 1.14 11.31 17.22 1.07 53 11
both log regression
log y=b0+b1log x x,y>0 y=exp b0 x**b1
ex8-2.stan
both log regression
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector<lower=0>[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~lognormal(b0+b1*log(x),s);
}
generated quantities{
vector[N] m0;
vector<lower=0>[N] y0;
for(i in 1:N){
m0[i]=b0+b1*log(x[i]);
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector<lower=0>[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*log(x1[i]);
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
qplot(log(x),log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -73.2126
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 13 -29.1697 0.00290829 0.000218058 1 1 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -29.17
## b0 1.28
## b1 1.89
## s 1.04
## m0[1] 3.89
## m0[2] 4.63
## m0[3] 2.94
## m0[4] 5.20
## m0[5] 5.60
## m0[6] 5.61
## m0[7] 4.05
## m0[8] 1.01
## m0[9] 4.59
## m0[10] -2.74
## m0[11] 2.77
## m0[12] 3.99
## m0[13] 5.48
## m0[14] 4.90
## m0[15] 5.26
## m0[16] 4.64
## m0[17] 5.14
## m0[18] 5.31
## m0[19] 5.14
## m0[20] 5.48
## y0[1] 39.80
## y0[2] 18.29
## y0[3] 9.75
## y0[4] 217.72
## y0[5] 63.03
## y0[6] 509.49
## y0[7] 85.65
## y0[8] 5.13
## y0[9] 53.25
## y0[10] 0.07
## y0[11] 42.10
## y0[12] 18.52
## y0[13] 87.21
## y0[14] 36.32
## y0[15] 124.60
## y0[16] 111.94
## y0[17] 135.48
## y0[18] 465.83
## y0[19] 53.36
## y0[20] 17.81
## m1[1] -2.74
## m1[2] 1.64
## m1[3] 2.85
## m1[4] 3.58
## m1[5] 4.11
## m1[6] 4.52
## m1[7] 4.86
## m1[8] 5.14
## m1[9] 5.39
## m1[10] 5.61
## y1[1] 0.01
## y1[2] 7.92
## y1[3] 14.91
## y1[4] 142.32
## y1[5] 157.01
## y1[6] 172.96
## y1[7] 157.79
## y1[8] 62.92
## y1[9] 260.96
## y1[10] 307.46mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -30.75 -30.40 1.32 1.09 -33.32 -29.31 1.00 891 933
## b0 1.26 1.27 0.46 0.45 0.50 2.00 1.00 488 595
## b1 1.90 1.90 0.26 0.25 1.49 2.35 1.00 610 738
## s 1.18 1.15 0.21 0.20 0.89 1.56 1.00 1299 1331
## m0[1] 3.89 3.89 0.27 0.27 3.44 4.32 1.00 941 1184
## m0[2] 4.63 4.64 0.28 0.28 4.15 5.09 1.00 1536 1275
## m0[3] 2.93 2.94 0.31 0.30 2.43 3.43 1.00 597 837
## m0[4] 5.21 5.22 0.32 0.30 4.67 5.73 1.01 2031 1355
## m0[5] 5.62 5.62 0.35 0.33 5.03 6.19 1.00 1749 1214
## m0[6] 5.62 5.63 0.35 0.33 5.04 6.20 1.00 1746 1214
## m0[7] 4.05 4.05 0.27 0.27 3.59 4.49 1.00 1049 1078
## m0[8] 0.99 1.00 0.49 0.48 0.17 1.77 1.00 486 595
## m0[9] 4.59 4.59 0.28 0.27 4.12 5.05 1.00 1497 1226
## m0[10] -2.78 -2.75 0.97 0.93 -4.44 -1.22 1.00 513 609
## m0[11] 2.76 2.76 0.32 0.31 2.23 3.27 1.00 569 764
## m0[12] 3.99 3.99 0.27 0.27 3.53 4.42 1.00 1006 1126
## m0[13] 5.49 5.49 0.34 0.32 4.93 6.05 1.00 1853 1250
## m0[14] 4.91 4.91 0.30 0.29 4.40 5.40 1.01 1861 1317
## m0[15] 5.27 5.27 0.32 0.31 4.71 5.80 1.00 2008 1355
## m0[16] 4.65 4.65 0.28 0.28 4.17 5.11 1.00 1550 1388
## m0[17] 5.15 5.16 0.31 0.30 4.62 5.67 1.01 2026 1375
## m0[18] 5.32 5.32 0.33 0.31 4.76 5.86 1.00 1978 1356
## m0[19] 5.15 5.16 0.31 0.30 4.62 5.66 1.01 2027 1375
## m0[20] 5.49 5.50 0.34 0.32 4.93 6.05 1.00 1849 1250
## y0[1] 108.47 46.43 345.78 46.27 6.19 350.67 1.00 1979 1958
## y0[2] 229.84 107.34 454.42 109.12 12.80 819.31 1.00 1982 1877
## y0[3] 42.76 18.47 113.03 18.91 2.38 140.63 1.00 1991 1932
## y0[4] 429.60 184.63 987.56 188.28 25.65 1466.08 1.00 1789 1737
## y0[5] 603.81 273.66 1149.40 276.96 32.56 2119.94 1.00 1886 1506
## y0[6] 600.09 269.02 1206.65 272.69 38.28 2078.40 1.00 2196 1873
## y0[7] 113.56 56.22 196.14 55.03 7.77 406.47 1.00 1832 2014
## y0[8] 6.98 2.83 20.30 2.98 0.34 22.00 1.00 1319 1670
## y0[9] 225.14 100.24 513.38 99.04 14.03 758.47 1.00 1882 1820
## y0[10] 0.25 0.06 1.37 0.08 0.01 0.80 1.00 897 1299
## y0[11] 34.37 15.62 71.60 15.72 2.06 110.92 1.00 1794 1908
## y0[12] 113.78 50.97 228.42 52.88 6.44 392.13 1.00 1783 1839
## y0[13] 631.91 232.53 2658.28 233.96 31.53 1979.47 1.00 2015 1866
## y0[14] 276.54 134.94 491.43 134.66 17.43 922.66 1.00 1981 1814
## y0[15] 438.80 189.07 1341.15 191.20 27.30 1389.58 1.00 1683 1786
## y0[16] 249.70 102.40 706.75 106.67 13.23 823.03 1.00 2027 1907
## y0[17] 355.59 170.31 632.65 171.38 22.80 1226.86 1.00 2205 1952
## y0[18] 444.03 197.32 997.61 201.42 26.88 1564.43 1.00 2146 1870
## y0[19] 346.25 159.91 611.51 162.48 21.89 1276.36 1.00 2270 1974
## y0[20] 543.73 232.44 1930.13 231.14 33.79 1597.37 1.00 1993 1795
## m1[1] -2.78 -2.75 0.97 0.93 -4.44 -1.22 1.00 513 609
## m1[2] 1.62 1.63 0.42 0.42 0.92 2.31 1.00 493 635
## m1[3] 2.84 2.85 0.31 0.31 2.33 3.35 1.00 582 832
## m1[4] 3.58 3.58 0.28 0.27 3.12 4.01 1.00 783 1141
## m1[5] 4.11 4.11 0.27 0.28 3.65 4.54 1.00 1091 1157
## m1[6] 4.52 4.53 0.28 0.28 4.05 4.98 1.00 1437 1243
## m1[7] 4.86 4.87 0.30 0.28 4.36 5.35 1.00 1803 1346
## m1[8] 5.15 5.16 0.31 0.30 4.62 5.67 1.01 2028 1375
## m1[9] 5.40 5.41 0.33 0.31 4.84 5.95 1.00 1925 1336
## m1[10] 5.62 5.63 0.35 0.33 5.04 6.20 1.00 1746 1214
## y1[1] 0.20 0.06 0.89 0.07 0.00 0.79 1.00 911 1193
## y1[2] 12.32 5.01 48.24 4.94 0.59 36.37 1.00 1644 1829
## y1[3] 39.36 17.96 95.58 17.90 2.36 125.22 1.00 1728 1722
## y1[4] 84.15 35.97 209.74 35.62 4.69 297.45 1.00 1924 1506
## y1[5] 123.15 60.15 207.89 59.94 7.91 426.91 1.00 1971 1856
## y1[6] 235.77 91.10 1683.98 91.17 12.89 703.81 1.00 2079 1886
## y1[7] 287.76 129.24 720.95 129.21 15.75 961.33 1.00 2034 1776
## y1[8] 385.82 181.41 863.84 186.17 22.73 1348.83 1.00 1966 1779
## y1[9] 510.92 223.51 1139.13 225.45 31.41 1705.22 1.00 1983 1703
## y1[10] 647.64 280.10 2096.45 282.79 35.76 2112.73 1.00 1875 1856
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
qplot(log(x),log(y),col=I('red'))+
geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=log(x),y=m),xy,col='black')
exponential increasing/decreasing
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)
b1>0 x->Infinity,y->Infinity
b1<0 x->Infinity,y->+0 n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)ex8-3-1.stan
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0*exp(b1*x),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*exp(b1*x[i]);
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*exp(b1*x1[i]);
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)## Initial log joint probability = -29.9215
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 15 7.95156 3.54027e-05 0.00348442 0.6984 0.6984 29
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ 7.95
## b0 10.14
## b1 -1.93
## s 0.41
## m0[1] 0.15
## m0[2] 0.16
## m0[3] 1.95
## m0[4] 3.33
## m0[5] 0.03
## m0[6] 6.49
## m0[7] 8.86
## m0[8] 0.08
## m0[9] 0.03
## m0[10] 0.17
## m0[11] 0.01
## m0[12] 0.00
## m0[13] 0.02
## m0[14] 0.00
## m0[15] 0.01
## m0[16] 0.22
## m0[17] 0.44
## m0[18] 0.01
## m0[19] 1.50
## m0[20] 0.00
## y0[1] -0.41
## y0[2] 0.24
## y0[3] 1.58
## y0[4] 3.18
## y0[5] 0.34
## y0[6] 6.59
## y0[7] 8.59
## y0[8] 0.19
## y0[9] 0.01
## y0[10] 0.10
## y0[11] -0.21
## y0[12] 0.53
## y0[13] -0.04
## y0[14] -0.08
## y0[15] 0.20
## y0[16] 0.01
## y0[17] 0.11
## y0[18] -0.56
## y0[19] 2.07
## y0[20] 0.01
## m1[1] 8.86
## m1[2] 3.26
## m1[3] 1.20
## m1[4] 0.44
## m1[5] 0.16
## m1[6] 0.06
## m1[7] 0.02
## m1[8] 0.01
## m1[9] 0.00
## m1[10] 0.00
## y1[1] 9.42
## y1[2] 4.28
## y1[3] 0.88
## y1[4] -0.16
## y1[5] 0.57
## y1[6] -0.08
## y1[7] 0.23
## y1[8] 0.38
## y1[9] -0.88
## y1[10] -0.17
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 9.23 9.63 1.46 1.12 6.28 10.77 1.01 535 778
## b0 10.19 10.19 0.56 0.55 9.27 11.12 1.01 1367 835
## b1 -1.96 -1.96 0.19 0.18 -2.29 -1.66 1.00 1502 1079
## s 0.47 0.45 0.09 0.07 0.35 0.63 1.03 331 393
## m0[1] 0.15 0.14 0.06 0.05 0.07 0.25 1.00 1605 1134
## m0[2] 0.16 0.15 0.06 0.05 0.07 0.27 1.00 1606 1134
## m0[3] 1.92 1.92 0.26 0.24 1.50 2.36 1.00 1847 1212
## m0[4] 3.28 3.29 0.29 0.26 2.80 3.74 1.00 2010 1240
## m0[5] 0.03 0.03 0.02 0.01 0.01 0.06 1.00 1572 1138
## m0[6] 6.47 6.47 0.28 0.26 6.02 6.93 1.00 2202 1452
## m0[7] 8.88 8.88 0.42 0.42 8.18 9.57 1.01 1497 974
## m0[8] 0.08 0.08 0.04 0.03 0.03 0.15 1.00 1589 1135
## m0[9] 0.04 0.03 0.02 0.02 0.01 0.07 1.00 1575 1139
## m0[10] 0.17 0.16 0.06 0.06 0.08 0.28 1.00 1608 1134
## m0[11] 0.01 0.01 0.01 0.01 0.00 0.03 1.00 1561 1158
## m0[12] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1549 1104
## m0[13] 0.02 0.02 0.01 0.01 0.01 0.05 1.00 1568 1138
## m0[14] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1545 1104
## m0[15] 0.01 0.01 0.01 0.00 0.00 0.02 1.00 1556 1158
## m0[16] 0.22 0.21 0.08 0.07 0.11 0.36 1.00 1618 1113
## m0[17] 0.44 0.43 0.12 0.11 0.26 0.66 1.00 1652 1127
## m0[18] 0.01 0.00 0.00 0.00 0.00 0.01 1.00 1555 1104
## m0[19] 1.48 1.47 0.24 0.22 1.10 1.88 1.00 1781 1231
## m0[20] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1545 1104
## y0[1] 0.16 0.15 0.49 0.46 -0.62 0.98 1.00 1955 1857
## y0[2] 0.18 0.16 0.48 0.47 -0.60 0.95 1.00 2157 1888
## y0[3] 1.91 1.91 0.55 0.52 1.01 2.78 1.00 1869 1923
## y0[4] 3.28 3.29 0.55 0.51 2.38 4.19 1.00 1609 1493
## y0[5] 0.04 0.04 0.47 0.45 -0.74 0.83 1.00 1952 1541
## y0[6] 6.46 6.45 0.56 0.52 5.53 7.34 1.00 2212 1592
## y0[7] 8.86 8.87 0.62 0.58 7.84 9.85 1.01 1712 1510
## y0[8] 0.08 0.09 0.48 0.45 -0.72 0.91 1.00 1888 1615
## y0[9] 0.02 0.04 0.48 0.45 -0.73 0.78 1.00 2053 1605
## y0[10] 0.17 0.17 0.48 0.47 -0.60 0.98 1.00 2039 1873
## y0[11] 0.01 0.01 0.49 0.48 -0.80 0.80 1.00 1891 1739
## y0[12] 0.01 0.00 0.47 0.45 -0.75 0.79 1.00 1797 1704
## y0[13] 0.03 0.04 0.48 0.44 -0.75 0.81 1.00 1852 1511
## y0[14] 0.00 0.00 0.46 0.45 -0.78 0.74 1.00 1913 1935
## y0[15] 0.01 0.00 0.48 0.46 -0.75 0.81 1.00 2118 1806
## y0[16] 0.23 0.23 0.49 0.45 -0.57 1.01 1.00 1984 1679
## y0[17] 0.45 0.44 0.49 0.47 -0.37 1.25 1.00 1915 1908
## y0[18] 0.02 0.02 0.48 0.45 -0.74 0.80 1.00 2028 1803
## y0[19] 1.46 1.46 0.53 0.51 0.59 2.34 1.00 1888 1749
## y0[20] -0.01 -0.02 0.48 0.45 -0.81 0.77 1.00 2123 1853
## m1[1] 8.88 8.88 0.42 0.42 8.18 9.57 1.01 1497 974
## m1[2] 3.21 3.22 0.29 0.26 2.73 3.67 1.00 1997 1231
## m1[3] 1.17 1.17 0.21 0.20 0.84 1.55 1.00 1742 1246
## m1[4] 0.43 0.42 0.12 0.11 0.25 0.65 1.00 1652 1127
## m1[5] 0.16 0.15 0.06 0.05 0.08 0.27 1.00 1607 1134
## m1[6] 0.06 0.06 0.03 0.02 0.02 0.12 1.00 1585 1135
## m1[7] 0.02 0.02 0.01 0.01 0.01 0.05 1.00 1569 1138
## m1[8] 0.01 0.01 0.01 0.00 0.00 0.02 1.00 1558 1158
## m1[9] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1551 1104
## m1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1545 1104
## y1[1] 8.87 8.87 0.61 0.59 7.87 9.86 1.00 1803 1654
## y1[2] 3.21 3.21 0.55 0.52 2.28 4.11 1.00 2248 1656
## y1[3] 1.17 1.17 0.53 0.51 0.31 2.02 1.00 2003 1565
## y1[4] 0.45 0.44 0.49 0.47 -0.30 1.24 1.00 1970 1784
## y1[5] 0.17 0.16 0.46 0.46 -0.56 0.94 1.00 1934 1879
## y1[6] 0.04 0.03 0.48 0.47 -0.72 0.86 1.00 1651 1800
## y1[7] 0.02 0.01 0.47 0.43 -0.74 0.80 1.00 2175 1817
## y1[8] 0.00 0.01 0.47 0.45 -0.78 0.76 1.00 2146 1525
## y1[9] 0.00 -0.02 0.46 0.43 -0.74 0.77 1.00 2057 1933
## y1[10] 0.00 -0.01 0.48 0.46 -0.76 0.79 1.00 1853 1846
ex8-3-2.stan
log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)
data{
int N;
vector<lower=0>[N] x;
vector<lower=0>[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0> b0;
real<lower=-10,upper=10> b1;
real<lower=0> s;
}
model{
y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=log(b0)+b1*x[i];
y0[i]=lognormal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=log(b0)+b1*x1[i];
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
qplot(x,log(y)),
ncol=2)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -233.923
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 24 -13.9674 0.000121639 0.000125591 0.9884 0.9884 37
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.mle## variable estimate
## lp__ -13.97
## b0 6.63
## b1 -1.91
## s 0.49
## m0[1] -7.55
## m0[2] -4.26
## m0[3] -6.12
## m0[4] -3.96
## m0[5] -3.09
## m0[6] -6.58
## m0[7] -0.27
## m0[8] -5.69
## m0[9] -6.76
## m0[10] -0.91
## m0[11] -2.66
## m0[12] -6.13
## m0[13] -3.52
## m0[14] -1.36
## m0[15] -7.24
## m0[16] 1.31
## m0[17] -5.29
## m0[18] -1.67
## m0[19] -4.83
## m0[20] -3.66
## y0[1] 0.00
## y0[2] 0.01
## y0[3] 0.00
## y0[4] 0.03
## y0[5] 0.07
## y0[6] 0.00
## y0[7] 0.25
## y0[8] 0.00
## y0[9] 0.00
## y0[10] 0.44
## y0[11] 0.02
## y0[12] 0.00
## y0[13] 0.03
## y0[14] 0.16
## y0[15] 0.00
## y0[16] 3.13
## y0[17] 0.01
## y0[18] 0.15
## y0[19] 0.00
## y0[20] 0.03
## m1[1] 1.31
## m1[2] 0.33
## m1[3] -0.66
## m1[4] -1.64
## m1[5] -2.63
## m1[6] -3.61
## m1[7] -4.60
## m1[8] -5.58
## m1[9] -6.57
## m1[10] -7.55
## y1[1] 4.33
## y1[2] 4.23
## y1[3] 0.52
## y1[4] 0.35
## y1[5] 0.17
## y1[6] 0.02
## y1[7] 0.01
## y1[8] 0.01
## y1[9] 0.00
## y1[10] 0.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "s" "m0" "y0" "m1" "y1"seeMCMC(mcmc,'m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.86 -12.51 1.42 1.16 -15.69 -11.35 1.00 612 600
## b0 7.84 7.28 2.95 2.32 4.39 12.83 1.01 330 452
## b1 -1.94 -1.94 0.10 0.10 -2.11 -1.78 1.01 378 605
## s 0.56 0.54 0.11 0.10 0.41 0.76 1.00 685 554
## m0[1] -7.59 -7.58 0.23 0.22 -7.97 -7.25 1.00 771 964
## m0[2] -4.25 -4.25 0.12 0.12 -4.46 -4.05 1.00 1675 1206
## m0[3] -6.13 -6.13 0.17 0.16 -6.42 -5.87 1.00 1327 1178
## m0[4] -3.94 -3.94 0.12 0.12 -4.15 -3.75 1.00 1165 1363
## m0[5] -3.06 -3.06 0.13 0.12 -3.28 -2.85 1.00 587 1047
## m0[6] -6.60 -6.60 0.18 0.18 -6.91 -6.32 1.00 1053 1179
## m0[7] -0.19 -0.20 0.23 0.22 -0.56 0.20 1.01 335 395
## m0[8] -5.70 -5.70 0.15 0.15 -5.96 -5.46 1.00 1719 1210
## m0[9] -6.79 -6.78 0.19 0.19 -7.11 -6.49 1.00 979 1137
## m0[10] -0.85 -0.85 0.21 0.20 -1.17 -0.51 1.01 344 466
## m0[11] -2.62 -2.62 0.14 0.13 -2.85 -2.39 1.00 469 747
## m0[12] -6.15 -6.15 0.17 0.16 -6.44 -5.89 1.00 1312 1178
## m0[13] -3.50 -3.50 0.13 0.12 -3.70 -3.30 1.00 811 1254
## m0[14] -1.30 -1.30 0.19 0.18 -1.60 -0.99 1.01 354 503
## m0[15] -7.27 -7.27 0.21 0.21 -7.63 -6.95 1.00 838 1031
## m0[16] 1.41 1.40 0.31 0.30 0.94 1.91 1.01 330 396
## m0[17] -5.29 -5.29 0.14 0.13 -5.53 -5.07 1.00 1827 1328
## m0[18] -1.62 -1.62 0.17 0.16 -1.90 -1.33 1.01 367 520
## m0[19] -4.83 -4.82 0.13 0.12 -5.05 -4.62 1.00 1888 1290
## m0[20] -3.64 -3.64 0.12 0.12 -3.85 -3.45 1.00 903 1224
## y0[1] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1687 1848
## y0[2] 0.02 0.01 0.01 0.01 0.01 0.04 1.00 1835 1739
## y0[3] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 2031 1879
## y0[4] 0.02 0.02 0.02 0.01 0.01 0.05 1.00 1832 1871
## y0[5] 0.06 0.05 0.04 0.02 0.02 0.12 1.00 2057 1999
## y0[6] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1985 1786
## y0[7] 1.00 0.84 0.71 0.49 0.31 2.19 1.00 1175 1568
## y0[8] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 2046 1854
## y0[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1708 1809
## y0[10] 0.51 0.42 0.40 0.23 0.16 1.17 1.01 1293 1425
## y0[11] 0.09 0.07 0.06 0.04 0.03 0.19 1.00 1943 1524
## y0[12] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1860 1823
## y0[13] 0.04 0.03 0.02 0.02 0.01 0.08 1.00 1899 1900
## y0[14] 0.33 0.27 0.24 0.15 0.10 0.73 1.00 1729 1624
## y0[15] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1539 1768
## y0[16] 5.12 3.98 4.64 2.35 1.40 11.88 1.00 860 1257
## y0[17] 0.01 0.01 0.01 0.00 0.00 0.01 1.00 1941 1746
## y0[18] 0.24 0.19 0.18 0.11 0.08 0.54 1.00 1814 1778
## y0[19] 0.01 0.01 0.01 0.00 0.00 0.02 1.00 1889 1579
## y0[20] 0.03 0.03 0.02 0.01 0.01 0.07 1.00 1959 1778
## m1[1] 1.41 1.40 0.31 0.30 0.94 1.91 1.01 330 396
## m1[2] 0.41 0.40 0.26 0.25 0.00 0.84 1.01 332 390
## m1[3] -0.59 -0.60 0.22 0.21 -0.93 -0.23 1.01 339 437
## m1[4] -1.59 -1.59 0.18 0.17 -1.87 -1.30 1.01 365 520
## m1[5] -2.59 -2.59 0.14 0.13 -2.82 -2.36 1.00 464 712
## m1[6] -3.59 -3.59 0.13 0.12 -3.80 -3.39 1.00 868 1188
## m1[7] -4.59 -4.59 0.13 0.12 -4.81 -4.38 1.00 1869 1238
## m1[8] -5.59 -5.59 0.15 0.15 -5.84 -5.35 1.00 1753 1258
## m1[9] -6.59 -6.59 0.18 0.18 -6.90 -6.30 1.00 1059 1179
## m1[10] -7.59 -7.58 0.23 0.22 -7.97 -7.25 1.00 771 964
## y1[1] 5.30 4.18 8.12 2.50 1.49 12.22 1.00 916 1557
## y1[2] 1.87 1.52 1.46 0.90 0.58 4.27 1.00 929 1085
## y1[3] 0.67 0.56 0.50 0.32 0.20 1.46 1.00 1398 1893
## y1[4] 0.25 0.21 0.17 0.11 0.08 0.54 1.00 1585 1551
## y1[5] 0.09 0.08 0.06 0.04 0.03 0.20 1.00 1863 1951
## y1[6] 0.03 0.03 0.02 0.01 0.01 0.07 1.00 1618 1508
## y1[7] 0.01 0.01 0.01 0.01 0.00 0.02 1.00 2166 1862
## y1[8] 0.00 0.00 0.00 0.00 0.00 0.01 1.00 1862 1796
## y1[9] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1861 1691
## y1[10] 0.00 0.00 0.00 0.00 0.00 0.00 1.00 1702 1703
y0=mcmc$draws('y0')
smy0=summary(y0)
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,log(y),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
growth curve
y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)ex8-4.stan
growth curve
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=10> b1;
real<lower=0> s;
}
model{
y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(1-exp(-b1*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(1-exp(-b1*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)## Initial log joint probability = -123.997
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -12.224 0.000253145 9.03431e-05 1 1 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -12.22
## b0 10.23
## b1 0.51
## s 1.12
## m0[1] 5.93
## m0[2] 7.06
## m0[3] 9.87
## m0[4] 7.04
## m0[5] 3.32
## m0[6] 9.02
## m0[7] 2.24
## m0[8] 9.44
## m0[9] 9.88
## m0[10] 9.19
## m0[11] 8.77
## m0[12] 10.05
## m0[13] 9.97
## m0[14] 9.99
## m0[15] 9.28
## m0[16] 8.72
## m0[17] 9.78
## m0[18] 9.81
## m0[19] 5.90
## m0[20] 9.96
## y0[1] 6.95
## y0[2] 7.06
## y0[3] 8.62
## y0[4] 8.13
## y0[5] 3.03
## y0[6] 9.95
## y0[7] 1.84
## y0[8] 9.71
## y0[9] 10.98
## y0[10] 8.56
## y0[11] 8.38
## y0[12] 10.80
## y0[13] 11.01
## y0[14] 8.08
## y0[15] 8.10
## y0[16] 9.14
## y0[17] 10.41
## y0[18] 9.58
## y0[19] 7.14
## y0[20] 11.82
## m1[1] 2.24
## m1[2] 5.01
## m1[3] 6.82
## m1[4] 8.00
## m1[5] 8.77
## m1[6] 9.28
## m1[7] 9.61
## m1[8] 9.82
## m1[9] 9.96
## m1[10] 10.05
## y1[1] 2.15
## y1[2] 4.81
## y1[3] 5.44
## y1[4] 8.88
## y1[5] 9.98
## y1[6] 9.92
## y1[7] 10.27
## y1[8] 7.91
## y1[9] 11.01
## y1[10] 12.02
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.22 -11.89 1.34 1.06 -14.83 -10.80 1.01 632 847
## b0 10.24 10.18 0.63 0.56 9.32 11.35 1.01 936 736
## b1 0.54 0.53 0.12 0.10 0.36 0.75 1.01 742 624
## s 1.27 1.24 0.23 0.21 0.95 1.70 1.00 845 929
## m0[1] 6.01 6.01 0.56 0.53 5.10 6.95 1.01 784 708
## m0[2] 7.11 7.12 0.51 0.49 6.24 7.93 1.01 825 768
## m0[3] 9.84 9.85 0.40 0.39 9.19 10.48 1.00 1478 1335
## m0[4] 7.10 7.10 0.52 0.49 6.23 7.92 1.01 824 768
## m0[5] 3.41 3.38 0.47 0.43 2.70 4.23 1.01 755 622
## m0[6] 9.01 9.01 0.33 0.32 8.46 9.54 1.00 1519 1161
## m0[7] 2.31 2.28 0.35 0.32 1.79 2.94 1.01 750 640
## m0[8] 9.41 9.42 0.32 0.32 8.87 9.93 1.00 2090 1357
## m0[9] 9.84 9.85 0.40 0.39 9.19 10.49 1.00 1466 1279
## m0[10] 9.17 9.17 0.32 0.31 8.63 9.69 1.00 1806 1422
## m0[11] 8.77 8.78 0.35 0.34 8.19 9.32 1.00 1278 1009
## m0[12] 10.02 10.02 0.47 0.45 9.27 10.83 1.00 1181 1081
## m0[13] 9.94 9.95 0.44 0.42 9.25 10.67 1.00 1310 1213
## m0[14] 9.95 9.96 0.44 0.42 9.25 10.69 1.00 1288 1149
## m0[15] 9.26 9.27 0.32 0.31 8.73 9.78 1.00 2046 1468
## m0[16] 8.72 8.73 0.35 0.35 8.13 9.27 1.00 1245 929
## m0[17] 9.74 9.75 0.37 0.36 9.13 10.34 1.00 1649 1463
## m0[18] 9.78 9.79 0.38 0.37 9.16 10.39 1.00 1586 1507
## m0[19] 5.99 5.98 0.56 0.53 5.08 6.92 1.01 784 708
## m0[20] 9.93 9.94 0.43 0.41 9.24 10.64 1.00 1331 1213
## y0[1] 6.08 6.06 1.43 1.38 3.82 8.50 1.00 1847 1833
## y0[2] 7.12 7.13 1.41 1.35 4.80 9.31 1.00 1693 1797
## y0[3] 9.81 9.83 1.35 1.28 7.64 12.04 1.00 1879 1759
## y0[4] 7.05 7.08 1.42 1.38 4.75 9.39 1.00 1584 1614
## y0[5] 3.43 3.45 1.38 1.32 1.18 5.63 1.00 1879 1721
## y0[6] 9.00 9.00 1.34 1.29 6.85 11.16 1.00 1846 1736
## y0[7] 2.32 2.30 1.33 1.22 0.15 4.57 1.00 1804 1502
## y0[8] 9.35 9.35 1.31 1.24 7.23 11.50 1.00 2082 1931
## y0[9] 9.88 9.86 1.34 1.25 7.79 12.04 1.00 1828 1775
## y0[10] 9.18 9.13 1.33 1.29 7.01 11.37 1.00 2305 1914
## y0[11] 8.78 8.74 1.35 1.27 6.55 11.10 1.00 2016 1914
## y0[12] 10.06 10.04 1.37 1.38 7.82 12.31 1.00 1628 1858
## y0[13] 9.93 9.95 1.36 1.35 7.68 12.08 1.00 1883 1893
## y0[14] 9.93 9.92 1.37 1.29 7.65 12.19 1.00 1958 1627
## y0[15] 9.30 9.30 1.34 1.26 7.10 11.51 1.00 2156 1566
## y0[16] 8.69 8.69 1.37 1.30 6.46 10.94 1.00 1837 1985
## y0[17] 9.78 9.80 1.32 1.28 7.57 11.93 1.00 1731 1859
## y0[18] 9.80 9.76 1.33 1.30 7.69 12.00 1.00 1956 2057
## y0[19] 6.01 6.00 1.42 1.37 3.71 8.25 1.00 1554 1798
## y0[20] 9.91 9.89 1.41 1.35 7.73 12.25 1.00 1849 1799
## m1[1] 2.31 2.28 0.35 0.32 1.79 2.94 1.01 750 640
## m1[2] 5.11 5.08 0.56 0.52 4.22 6.07 1.01 769 644
## m1[3] 6.88 6.89 0.53 0.50 6.00 7.74 1.01 812 731
## m1[4] 8.03 8.04 0.43 0.41 7.29 8.69 1.01 939 827
## m1[5] 8.77 8.78 0.35 0.34 8.19 9.32 1.00 1278 1009
## m1[6] 9.25 9.26 0.32 0.31 8.72 9.78 1.00 2038 1491
## m1[7] 9.57 9.58 0.34 0.33 9.01 10.11 1.00 1933 1355
## m1[8] 9.79 9.80 0.38 0.37 9.16 10.41 1.00 1567 1484
## m1[9] 9.93 9.94 0.43 0.41 9.24 10.64 1.00 1329 1213
## m1[10] 10.02 10.02 0.47 0.45 9.27 10.83 1.00 1181 1081
## y1[1] 2.28 2.27 1.34 1.30 0.17 4.52 1.00 1892 1767
## y1[2] 5.09 5.04 1.41 1.31 2.82 7.46 1.00 1610 1775
## y1[3] 6.91 6.91 1.39 1.39 4.71 9.24 1.00 1617 1600
## y1[4] 8.03 8.05 1.38 1.27 5.73 10.29 1.00 1791 1700
## y1[5] 8.78 8.78 1.30 1.26 6.62 10.86 1.00 1953 1654
## y1[6] 9.24 9.27 1.31 1.25 7.05 11.36 1.00 2051 1944
## y1[7] 9.57 9.58 1.31 1.28 7.40 11.71 1.00 1876 1884
## y1[8] 9.79 9.78 1.37 1.30 7.56 12.09 1.00 2093 1805
## y1[9] 9.94 9.95 1.36 1.32 7.69 12.13 1.00 1700 1811
## y1[10] 10.01 9.98 1.33 1.30 7.84 12.25 1.00 1664 1753
sigmoid curve
y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)ex8-5.stan
sigmoid curve
data{
int N;
vector[N] x;
int Ym;
array[N] int y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=100> b0;
real<lower=0,upper=100> b1;
}
model{
y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
array[N] int y0;
for(i in 1:N){
y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
}
array[N1] int y1;
for(i in 1:N1){
y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
}
}n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-5.stan')
mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc$metadata()$stan_variables## [1] "lp__" "b0" "b1" "y0" "y1"seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -44.74 -44.44 1.01 0.71 -46.81 -43.80 1.00 456 914
## b0 3.86 3.86 0.12 0.12 3.67 4.07 1.00 1195 1004
## b1 2.39 2.33 0.49 0.45 1.70 3.33 1.01 352 340
## y0[1] 9.96 10.00 0.21 0.00 10.00 10.00 1.00 1655 NA
## y0[2] 4.46 4.00 1.70 1.48 2.00 7.00 1.00 2055 1935
## y0[3] 9.67 10.00 0.63 0.00 8.00 10.00 1.00 1240 NA
## y0[4] 10.00 10.00 0.00 0.00 10.00 10.00 NA NA NA
## y0[5] 10.00 10.00 0.05 0.00 10.00 10.00 1.00 2027 NA
## y0[6] 5.38 5.00 1.74 1.48 3.00 8.00 1.01 1430 1911
## y0[7] 7.40 8.00 1.58 1.48 5.00 10.00 1.00 1337 NA
## y0[8] 3.05 3.00 1.54 1.48 1.00 6.00 1.00 1842 1860
## y0[9] 2.31 2.00 1.41 1.48 0.00 5.00 1.00 1853 1902
## y0[10] 0.01 0.00 0.12 0.00 0.00 0.00 1.00 1834 1834
## y0[11] 0.44 0.00 0.66 0.00 0.00 2.00 1.00 1881 1915
## y0[12] 10.00 10.00 0.06 0.00 10.00 10.00 1.00 2025 NA
## y0[13] 10.00 10.00 0.03 0.00 10.00 10.00 1.00 2018 NA
## y0[14] 10.00 10.00 0.04 0.00 10.00 10.00 1.00 2019 NA
## y0[15] 1.74 2.00 1.28 1.48 0.00 4.00 1.00 1505 1841
## y0[16] 0.04 0.00 0.21 0.00 0.00 0.00 1.00 1803 1774
## y0[17] 0.01 0.00 0.10 0.00 0.00 0.00 1.00 2057 2057
## y0[18] 2.28 2.00 1.42 1.48 0.00 5.00 1.00 2010 1970
## y0[19] 0.31 0.00 0.58 0.00 0.00 1.00 1.00 1712 1693
## y0[20] 0.03 0.00 0.17 0.00 0.00 0.00 1.00 1931 1954
## y1[1] 0.01 0.00 0.11 0.00 0.00 0.00 1.00 2065 2065
## y1[2] 0.07 0.00 0.27 0.00 0.00 1.00 1.00 1954 1972
## y1[3] 0.41 0.00 0.66 0.00 0.00 2.00 1.00 1389 1413
## y1[4] 2.58 2.00 1.46 1.48 0.00 5.00 1.00 1739 1962
## y1[5] 7.39 8.00 1.57 1.48 5.00 10.00 1.00 1404 NA
## y1[6] 9.52 10.00 0.77 0.00 8.00 10.00 1.00 1118 NA
## y1[7] 9.92 10.00 0.30 0.00 9.00 10.00 1.00 1748 NA
## y1[8] 9.98 10.00 0.12 0.00 10.00 10.00 1.00 2080 NA
## y1[9] 9.99 10.00 0.07 0.00 10.00 10.00 1.00 1267 NA
## y1[10] 10.00 10.00 0.04 0.00 10.00 10.00 1.00 2020 NA
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=ymed),xy,col='black')
up and down
y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))ex8-6.stan
up and down
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real<lower=0,upper=200> b0;
real<lower=0,upper=1> b1;
real<lower=0,upper=1> b2;
real<lower=0,upper=10> s;
}
model{
y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
y1[i]=normal_rng(m1[i],s);
}
}n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-6.stan')
fn(mdl,data)## Initial log joint probability = -5547.9
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## 54 -7.43127 1.90425e-05 0.000920521 1 1 118
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -7.43
## b0 95.89
## b1 0.03
## b2 0.21
## s 0.88
## m0[1] 56.55
## m0[2] 51.35
## m0[3] 45.09
## m0[4] 25.45
## m0[5] 60.17
## m0[6] 48.96
## m0[7] 30.89
## m0[8] 30.35
## m0[9] 60.43
## m0[10] 54.36
## m0[11] 30.60
## m0[12] 58.95
## m0[13] 59.20
## m0[14] 53.62
## m0[15] 40.88
## m0[16] 36.12
## m0[17] 60.56
## m0[18] 27.90
## m0[19] 49.48
## m0[20] 26.63
## y0[1] 54.94
## y0[2] 50.94
## y0[3] 46.29
## y0[4] 26.44
## y0[5] 59.58
## y0[6] 48.06
## y0[7] 30.93
## y0[8] 30.34
## y0[9] 60.29
## y0[10] 53.52
## y0[11] 31.04
## y0[12] 58.90
## y0[13] 58.66
## y0[14] 53.37
## y0[15] 43.07
## y0[16] 34.99
## y0[17] 60.97
## y0[18] 26.25
## y0[19] 50.06
## y0[20] 26.61
## m1[1] 54.36
## m1[2] 60.57
## m1[3] 58.27
## m1[4] 53.23
## m1[5] 47.60
## m1[6] 42.19
## m1[7] 37.25
## m1[8] 32.83
## m1[9] 28.91
## m1[10] 25.45
## y1[1] 54.67
## y1[2] 60.44
## y1[3] 59.99
## y1[4] 53.23
## y1[5] 46.92
## y1[6] 42.42
## y1[7] 36.12
## y1[8] 32.95
## y1[9] 28.47
## y1[10] 23.91
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.3 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.31 -11.03 1.43 1.36 -13.98 -9.54 1.01 472 813
## b0 95.93 95.80 2.76 2.53 91.68 100.52 1.01 495 738
## b1 0.03 0.03 0.00 0.00 0.03 0.03 1.01 588 838
## b2 0.21 0.21 0.01 0.01 0.20 0.23 1.01 525 773
## s 1.01 1.00 0.19 0.18 0.70 1.35 1.15 21 28
## m0[1] 56.54 56.53 0.37 0.37 55.92 57.13 1.02 720 971
## m0[2] 51.34 51.35 0.36 0.35 50.76 51.93 1.01 680 1006
## m0[3] 45.08 45.09 0.30 0.30 44.58 45.57 1.01 1546 1293
## m0[4] 25.46 25.46 0.41 0.40 24.78 26.14 1.01 897 1365
## m0[5] 60.15 60.15 0.43 0.40 59.45 60.88 1.02 1706 1646
## m0[6] 48.96 48.96 0.34 0.33 48.40 49.51 1.01 816 1185
## m0[7] 30.90 30.90 0.34 0.33 30.33 31.46 1.01 1291 1564
## m0[8] 30.36 30.36 0.35 0.34 29.78 30.93 1.01 1235 1600
## m0[9] 60.41 60.41 0.35 0.34 59.84 60.98 1.02 1885 1322
## m0[10] 54.37 54.36 0.71 0.67 53.20 55.54 1.01 842 1342
## m0[11] 30.61 30.62 0.35 0.34 30.04 31.18 1.01 1260 1600
## m0[12] 58.93 58.92 0.36 0.36 58.34 59.51 1.02 1114 1045
## m0[13] 59.20 59.19 0.50 0.49 58.36 60.04 1.02 1443 1386
## m0[14] 53.60 53.60 0.37 0.37 52.98 54.21 1.02 659 1054
## m0[15] 40.88 40.87 0.28 0.27 40.42 41.33 1.01 1997 1639
## m0[16] 36.12 36.13 0.29 0.29 35.65 36.60 1.01 1833 1846
## m0[17] 60.55 60.54 0.36 0.35 59.96 61.13 1.01 1972 1395
## m0[18] 27.91 27.91 0.38 0.36 27.28 28.53 1.01 987 1576
## m0[19] 49.47 49.47 0.34 0.34 48.91 50.02 1.01 771 1162
## m0[20] 26.64 26.64 0.40 0.38 25.98 27.29 1.01 938 1368
## y0[1] 56.52 56.53 1.13 1.09 54.68 58.30 1.01 1868 1718
## y0[2] 51.32 51.32 1.11 1.08 49.50 53.09 1.00 1762 1874
## y0[3] 45.08 45.08 1.06 0.99 43.32 46.84 1.01 1934 1805
## y0[4] 25.43 25.43 1.12 1.06 23.59 27.22 1.00 1676 1370
## y0[5] 60.12 60.15 1.09 1.04 58.29 61.86 1.00 1890 1798
## y0[6] 48.95 48.97 1.14 1.08 47.05 50.79 1.00 2008 1749
## y0[7] 30.89 30.86 1.09 1.04 29.11 32.70 1.01 2033 1731
## y0[8] 30.36 30.39 1.09 1.07 28.56 32.08 1.01 1889 1775
## y0[9] 60.44 60.43 1.06 1.00 58.78 62.17 1.01 1598 1619
## y0[10] 54.35 54.31 1.25 1.20 52.29 56.40 1.00 1514 1543
## y0[11] 30.63 30.63 1.13 1.10 28.84 32.46 1.00 1828 1769
## y0[12] 58.97 58.96 1.11 1.11 57.14 60.78 1.01 2078 1585
## y0[13] 59.18 59.19 1.15 1.07 57.21 61.11 1.01 1768 1663
## y0[14] 53.65 53.64 1.12 1.05 51.80 55.54 1.00 2040 1696
## y0[15] 40.87 40.90 1.06 1.03 39.11 42.59 1.00 1693 1585
## y0[16] 36.12 36.14 1.07 1.04 34.38 37.84 1.01 1848 1662
## y0[17] 60.57 60.55 1.08 1.07 58.88 62.40 1.01 1708 1865
## y0[18] 27.91 27.89 1.09 1.06 26.13 29.69 1.00 1801 1808
## y0[19] 49.46 49.44 1.08 1.05 47.68 51.35 1.01 1932 1823
## y0[20] 26.65 26.67 1.10 1.04 24.81 28.40 1.00 1896 1677
## m1[1] 54.37 54.36 0.71 0.67 53.20 55.54 1.01 842 1342
## m1[2] 60.56 60.56 0.38 0.36 59.93 61.18 1.01 1960 1533
## m1[3] 58.26 58.24 0.36 0.36 57.66 58.84 1.02 929 1004
## m1[4] 53.21 53.21 0.37 0.36 52.60 53.82 1.02 658 923
## m1[5] 47.59 47.59 0.33 0.32 47.06 48.13 1.01 963 1211
## m1[6] 42.19 42.19 0.28 0.28 41.71 42.65 1.01 1882 1497
## m1[7] 37.25 37.25 0.28 0.28 36.78 37.71 1.02 1942 1809
## m1[8] 32.83 32.83 0.32 0.31 32.30 33.37 1.01 1472 1658
## m1[9] 28.92 28.92 0.37 0.36 28.31 29.52 1.01 1086 1571
## m1[10] 25.46 25.46 0.41 0.40 24.78 26.14 1.01 897 1365
## y1[1] 54.37 54.37 1.20 1.16 52.43 56.30 1.01 1676 1704
## y1[2] 60.58 60.58 1.09 1.04 58.82 62.36 1.01 2040 1825
## y1[3] 58.26 58.26 1.10 1.07 56.47 60.05 1.01 1955 1357
## y1[4] 53.21 53.20 1.05 1.01 51.53 54.94 1.01 2089 1792
## y1[5] 47.60 47.59 1.05 1.00 45.87 49.29 1.00 1895 1683
## y1[6] 42.17 42.16 1.06 1.00 40.50 43.90 1.01 2036 1750
## y1[7] 37.26 37.26 1.12 1.08 35.40 39.10 1.00 1977 1361
## y1[8] 32.80 32.81 1.08 1.02 31.03 34.51 1.00 1928 1768
## y1[9] 28.91 28.92 1.12 1.09 27.03 30.64 1.01 1773 1678
## y1[10] 25.46 25.46 1.09 1.03 23.61 27.29 1.01 1537 1648
